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Abstract 

A simple model of the driven motion of interacting particles in a two di- 
mensional random medium is analyzed, focusing on the critical behavior near 
to the threshold that separates a static phase from a flowing phase with a 
steady-state current. The critical behavior is found to be surprisingly robust, 
being independent of whether the driving force is increased suddenly or adia- 
batically. Just above threshold, the flow is concentrated on a sparse network 
of channels, but the time scale for convergence to this fixed network diverges 
with a larger exponent that that for convergence of the current density to its 
steady-state value. This is argued to be caused by the "dangerous irrelevance" 
of dynamic particle collisions at the critical point. Possible applications to 
vortex motion near to the critical current in dirty thin film superconductors 
are discussed briefly. 

PACS number(s): 05.60.+W, 74.60.Ge, 62.20.Fe 
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I. INTRODUCTION 



Vortices in superconductors that are driven by an applied electric current and pinned by 
disorder exhibit a wide range of interesting behavior. These are examples of non-linear trans- 
port in random media which is collective in that the interactions between the transported 
can roughly be divided into two classes: "elastic" and "plastic" depending on whether or not 
the interactions between the transported particles are strong enough to maintain the parti- 
cles in an extended elastic structure as they move. Our interest here is plastic flow where the 
interactions between the particles and the random medium ("pinning") are strong enough 
to break up any elastic structure. We will focus on the very strong pinning limit for which 
"channel flow" can occur: where the flow is not only plastic but dominated by particles 
moving along a sparse network of persistent channels. 

Plastic flow of vortices has attracted a lot of recent interest .!'! A number of experimental 
measurements in 2H-NbSe2 have been attributed to plastic flow, including an unexpected 
dip in the electrical resistance just below H c2 ("peak effect"),!'! anomalous I-V curvesjl 
generation of unusual broadband noise,! small angle neutron scattering measurements^ and 
fingerprint phenomena where the detailed shape of the I-V curve is repeatable for a single 
sample but differs between samples.! Similar phenomena have recently been observed in 
two-dimensional amorphous Mo7yGe 2 3 films! which are supported by numerical simulations 
for two-dimensional systemdlUll which clearly see vortex motion dominated by flow along 
narrow "filamentary" channels. In addition to this indirect evidence, realtime images of 
moving vortices in thin films have been recorded.EHlli which clearly show individual vortices 
moving along narrow paths of least resistance. 

Some of these experiments also suggest that in the regimes studied plastic flow only 
persists for forces just above the threshold for the onset of motion, i.e., near the critical 
current. At higher current the vortex lattice appears to become more ordered. Previous 
theoretical wor k§E that studied plastic flow has primarily concentrated on the ordering 
and break up of an elastic lattice by the randomness, concluding that, at least for small 
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driving forces in two dimensions, the lattice will always break up and the flow will become 
plastic. Here we take a complementary approach that considers the extreme limit of a fully 
plastic state with only hard-core intervortex interactions. 

In this paper we extend work on a simple modeffi — roughly of point vortices in a thin 
film. It exhibits two phases: if the driving force is small the vortices are all trapped and 
there is no steady-state current, but if the force exceeds a finite threshold the vortices move 
in a static channel network whose configuration is determined by the pinning in the sample. 
In the steady-state just above threshold, the channels are far apart but each channel carries 
a high vortex current density. 

We focus on the dynamic critical behavior near threshold and on the development of the 
channel network above threshold. This study will be primarily numerical aided by scaling 
analysis. The critical behavior appears to represent a universality class of non-equilibrium 
"phase" transitions that probably also includes some aspects of a continuum fluid model 
studied by Narayan and Fisher.0 



A. Outline 



In the rest of the Introduction we briefly explain the model system we have investigated 
and then outline our main results. In Sec. || we provide a more detailed description of the 



model. The dynamic critical phenomena and scaling properties are explained in Sec. |TJ, 
while the development of the channel network is studied in Sec. [IV]. Finally we compare our 
results to related work and consider possible applications in Sec. [V]. 



B. Model 

The model studied here was introduced in Ref. [H] (where it was called the "one-deep" 
model). We briefly summarize its features here and define it in more detail in Sec. [TT]. 

This simple model for the motion of vortices in a two-dimensional random medium, 
treats the vortices as point particles with a hard-core interaction which move on a lattice. 
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The randomness is represented by each site's capacity to trap arriving particles and by the 
existence of a locally preferred direction — the "primary outlet" — for motion out of each site 
(see Fig. [I]). The local capacity and preferred direction can be thought of as barriers for a 
particle to leave by one of the site's outlets. Particles only leave through a particular outlet 
of a site when the number of particles present exceeds the barrier height. Increasing the 
driving force (F) corresponds to lowering all of the barrier heights which causes some sites 
to overflow thereby releasing particles. If one excess particle is present on a site, it always 
leaves thought the primary outlet. If more than one excess particle is present, one will leave 
through each outlet — a "split". Decreasing the force raises barriers, creating more traps and 
increasing the capacity of existing traps. 

Various possible histories of the system are represented by the way in which the force is 
changed, in this paper we will mainly consider increasing the force very quickly ("sudden" 
forcing) and then waiting for the particles to reach a steady-state. But we will also investigate 
what happens when the force is increased very slowly ( "adiabatic" forcing) for comparison. 

C. Dynamic critical behavior 

The behavior of the system is divided into two phases: A flowing phase where a steady- 
state current flows and a static phase with no current. In the static phase, the system 
has two types of sites: saturated sites at which any additional particle would overflow, and 
unsaturated sites or "traps". In the static phase, the particles will respond to a sudden 
increase in the force only in a transient manner, with some fraction of the particles mov- 
ing varying distances downhill before being trapped at previously unsaturated sites. The 
static configuration can be probed by adding a single particle, which, if it is added to a 
saturated site, will flow through primary outlets down a "drainage tree" until it stops at 
the unsaturated site at the bottom of the tree. At a critical force F c , the lengths of some of 
the trees will diverge and above the critical force the system will have steady-state flow. A 
fundamental property of our model is that at large times above threshold the particle flow 
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converges to a fixed channel network with each channel carrying its maximum capacity.El 
This implies that there are three types of lattice sites in the steady-state: network sites that 
carry moving particles, and off-network sites with no particle flow which can be either in 
the drainage basin (sites that drain to the channel network), or in finite drainage trees. The 
network sites form persistent channels that drain the system. The channels can join and 
split as they move through the system. In an infinite system, which sites are on the network 
(for a particular current) depends only on which outlets are the primary outlets and not on 
the initial particle placements. 

We now summarize the critical behavior. In the flowing phase in a large system the 
steady-state current density (i.e., the mean number of particles per lattice site), after a 
sudden increase in the force to a final value of F just above F c , is found to be 

J~(F-F c f. (1) 

Another quantity of interest in the flowing phase is the fraction of sites in the drainage 
basin, i.e., sites to which an added particle, following primary outlets, would continue to 
move indefinitely and become part of the steady-state current. The fraction of sites in the 
basin is found to increase as 

B~(F- F c ) r . (2) 

Typical results from a set of simulations showing this behavior are shown in Fig. 0. 
Numerical simulations of the model yield values of these exponents 

P = 1.53 ±0.03, (3) 
r = 0.49 ±0.02. (4) 

These error bars, and all uncertainties quoted and plotted here, are estimates of statistical 
l-a errors. More details of the accuracy and independence of these measurements are given 
in Sec. [Tj. 
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1. Length scales 



We are, of course, primarily interested in the behavior of infinite systems. But numerical 
studies of a range of system sizes yield both ways to extrapolate to infinite systems and 
information on the characteristic length scales in the system. Our scaling is anisotropic in 
the two directions, we will generally measure length (L) in the downhill direction and widths 
(W) in the horizontal direction. 

The behavior of a finite system of length L is controlled by a correlation length as 
described in Sec. [TT1]. This length scale controls the finite size scaling and is found to diverge 
with an exponent v as F — > F c . In order to obtain the critical force and the correlation 
length it is convenient to work with H(L, W, F) defined as the probability that there is a 
steady-state current for a system of length L and width W at force F. Using the scaling of 
IT we find 

v= 1.62 ±0.04 (5) 

on both sides of the threshold. 

Below threshold at long times all particles are at rest and the sites can be designated 
saturated and unsaturated ( "traps" ) depending on whether an added particle will flow out of 
the site. The saturated sites and their primary outlets (see Fig. ^ form connected drainage 
trees such that an excess particle added to a site on a drainage tree comes to rest on the 
unsaturated terminus site at the bottom of the tree. The only important length scale is the 
characteristic tree length above which trees are exponentially unlikely, this we identify as 
proportional to £j. 

However, above threshold the situation is rather more complicated: there are two length 
scales present in the channel network. The characteristic vertical distance scale of the 



channel network backbone, £ n , which, as shown in Ref. 21, diverges as 



& ~ jt (6) 
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as J — > 0; and the length which controls finite size scaling If £ n was identified with 

£f then we would expect v — 2/3 which is clearly ruled our by the data; instead we find 
(3/v — 0.96 ± 0.02. Nevertheless, £ n definitely controls various physical quantities which 
suggests that both lengths are important above threshold. 

The presence of two different correlation lengths above threshold is rather unusual (al- 
though it has appeared in other non-equilibrium dynamics critical points, for example sliding 
charge density waves.il) In our system the two lengths control two different equilibration 
processes. The decay of current transients is rapid with characteristic time £/, however redis- 
tribution of the current to the steady-state channel network is much slower. The convergence 
to a fixed channel network can be seen by studying the distribution of time-averaged local 
currents through different sites of the lattice. For large systems at long times most of the 
sites either always contain moving particles or never contain moving particles. This distri- 
bution has been studied numerically and is is indeed found to converge to two delta functions 
(at currents of zero and one) as the system size increases with characteristic time scale of 
order £ n and a power law tail (1/t 1 / 4 ) at long times. Because this approach to a steady-state 
pattern is so slow near threshold (and £ n ^> £/) experiments or simulations may not reach 
the steady-state of the current network even if the current itself does reach its steady-state 
value. 

The forcing in our model is strongly anisotropic and particles move much further parallel 
to the force (downhill) than perpendicular to it. The correlation lengths £ f and £ n are defined 
parallel to the force. For perpendicular (horizontal) distances we define two characteristic 
lengths £ f x and £ n± which we expect to behave as 

Ox ~ # (7) 

and 

finx~£-- (8) 



From Ref. [21], the properties of the channel network fix 



«n = 2" ( 9 ) 

Numerically, we obtain 

a = 0.50 ±0.05. (10) 



2. Adiabatic forcing and universality class 

Our model seems to represent a new universality class. Some of its characteristics are 
reminiscent of directed invasion percolation but the behavior is qualitatively different. An 
important feature of directed invasion percolationl! is that a local advance of "fluid" is not 
retarded by any elastic interaction with the rest of the fluid (unlike in models with elasticity 
representing surface tension). Our model goes further because the advancement of a local 
structure ahead of the rest of the particles is actively favored: as the bottom of a drainage 
tree moves down it collects more particles from above and hence is more likely to advance 
further under the "weight" of particles from above. 

The exponents z/, a and T have direct analogs in directed percolation, where they take 
the values^ 

^dir-perc ~ 1.73 (11) 
"dir-perc « 0.63 (12) 

/Wrc « 0-28 (13) 

with /^dir-perc analogous to our T. These are significantly different from our values. In- 
deed, even in mean field theory, which has only been analyzable for adiabatic forcing below 
threshold, our system is already different from directed percolation. 

An important question in investigating a new class of critical phenomena is how many 
independent exponents there are. We will argue that the scaling exponents of the current 
density and of the drainage basin should be related, 
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P = l + T. (14) 

If in fact this is an equality, as is supported by the numerics and analytic arguments, then 
there would be only three independent exponents v and a) as for directed percolation. 
In fact, we will argue that the exponent a should be exactly given by 

a = i, (15) 

(consistent with the numerics) which would imply only two independent exponents in our 
system. These issues are closely tied to the relationship between the behavior with adiabatic 
versus sudden forcing. 

With adiabatic forcing only one particle moves at a time when the system is below 
threshold. This means that there are no particle collisions — i.e., two particles never arrive 
at a site at the same time — and all particles move through primary outlets, i.e., there are 
no splits below threshold. In contrast, for sudden forcing many particles are released at 
once and then move down the lattice until they find traps. Even below threshold (where 
there is no steady-state current) there will be a transient flow with many collisions between 
particles. These collisions cause particles to use secondary outlets and explore sites which 
would otherwise not have been reached. This means extra traps are found and the threshold 



force is increased. Indeed we find (Sec. |III B|) critical forces for the two cases that are quite 
close but clearly different. But a more subtle question is whether the collisions change the 
universality class. 

If the only effect of the collisions is for a finite fraction of the transient moving particles 
to move along different outlets then the type of forcing should not effect the universality 
class, as argued in Sec. [TirT]. For adiabatic forcing the scaling relation 

(3 a = l + T a , (16) 

obtained by considering the effects of extra added particles, should be exact (Sec. |lllFj ). If 
the effect of the collisions is unimportant asymptotically then the two types of forcing are 
in the same universality class and the scaling relation [Eqs. [14|] holds. 



To test whether the universality classes are different one might hope to measure different 
critical exponents for the adiabatic case. Unfortunately direct simulations above threshold 
are difficult because, in principle, the system must be allowed to reach a steady-state before 
adding each new particle when above threshold. Thus we are restricted to scaling below 
threshold. From this we find 

v a = 1.60 ±0.05 (17) 

which overlaps with the value for the sudden case. We also find T a /v a = 0.29 ± 0.02 which 
when combined with the above value for v a leads to 

r a = 0.49 ± 0.04, (18) 

which again is well within the error bars of the sudden forcing value. 

On the basis of the data and the physical arguments of Sec. [II IF] we conjecture that 



P = 1 + T = 1.51 ±0.03 (19) 

is exact for both sudden and adiabatic forcing and that they are in the same universality 
class. 

This conclusion, and the presence of two diverging correlation lengths in the flowing 
phase with only the longer one manifesting the effects of collisions suggests that at the 
critical points, collisions between moving particles and the subsequent divergence of their 
paths are dangerously irrelevant. The collisions are clearly crucial at long times above 
threshold: if they did not occur at all the moving particles would eventually collect together 
on a single drainage tree with a diverging local current density. But near threshold the 
physics of the collisions appears to manifest itself only on very long length (or time) scales, 
£ n . Most of the physical properties near threshold — the mean current density, the drainage 
basin density, the statistics of the drainage trees etc — do not depend on the collisions. In 
Sec. [IV D| we will use the behavior in the collisionless regime above threshold on lengths 
scales, L, in the range 
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£/ < L < a, (20) 

to argue that the exponent a should be exactly 1/2. This will also make natural contact 
with the continuous fluid river model system studied earlier .0 



II. MODEL 

We now define the details of the model studied here. The model consists of particles 
moving on a lattice. The (quenched) randomness of the medium is represented by how many 
particle can be held (trapped) at a site before further particles "overflow" to the next site, 
and a local rule for distributing overflowing particles among nearest neighbor downhill sites. 
Our main focus is on two dimensional systems and we consider motion on a square lattice 
oriented as shown in Fig. [I]. The applied force acts along one of the diagonals and particles 
can move to nearest neighbor sites along the two possible downhill directions, so that each 
lattice site has two inlets and two outlets. 

The way in which particles move out of a site where the local capacity is exceeded 
depends on the exact variant of the model under consideration. In this paper we use only 
the simplest case (the "one-deep model") which is described in detail in the next section. 
However, conceptually it is useful to think more generally in terms of choosing a barrier 
height for each of the two outlets from some distribution. When the number of particles is 
less than the lower outlet barrier, then no particles leave. If the number of particles at the 
site exceeds only one of the outlet heights then all of the particles above the lower outlet's 
height leave by that outlet (which we call the "primary outlet"). If the number of particles 
exceeds the heights of both outlets then the number above the higher outlet is divided in 
some way between the two outlets. 

It is convenient to work with integer outlet barriers that are the least integer greater 
than the continuous valued barrier height. Increasing the driving force, F, corresponds to 
lowering all of barrier heights uniformly, resulting in some fraction of the integer barrier 
heights decreasing by one. This is equivalent to adding a particle to some fraction of the 
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sites. However, the particle sites are not chosen completely independently, since decreasing 
all barriers together means that a new particle will not appear a second time on any site 
until a new one has appeared on every site. For most of our purposes this difference will be 
unimportant. 

A. One deep model 

The simplest version of this model was introduced in Ref. [2l|, it allows at most one 
particle to past through each outlet in one time step (hence "one-deep"). The horizontal 
rows of the rotated square lattice of Fig. [I] are denoted by y = 1,2,3..., numbering from 
the top down, and the sites on row y by (x,y), with x + y an even integer. The (integer) 
number of particles n(x,y,t) on site (x, y) at time t is measured relative its capacity, i.e., 
to the lower integer outlet barrier so that n < implies that the site is unsaturated, i.e., 
another particle can be trapped while n = implies the site is saturated and n > that the 
site will overflow. 

If n(x, y, t) — 1 then at time t + 1, one particle is moved to the site in the next lower row 
to which the randomly chosen (but fixed) primary outlet of the site connects, i.e., to one of 
the sites {x ± 1, y + 1). If n(x, y, t) = 2 then at time t + 1, one particle is moved to each of 
the two sites (x±l,y + 1). These rules ensure that if we start with all n(x, y, t = 0) < 2 and 
update all sites simultaneously then no more than one particle passes through any outlet at 
any time step. The model is further simplified if we also restrict unsaturated sites (traps) 
to a capacity of one added particle so that 

-l<n(x,y,t)<2. (21) 

A variety of lattice sizes were chosen for numerical simulations. Lattice lengths (L) are 
measured on the y-scale of Fig. [I], widths (W) on the x-scale. A directed random walker 
on the lattice moving along outlets has a diffusion constant of 1/2, i.e., the mean square 
displacement of a walk of length y grows as (Ax(y) 2 ) ~ y for large y. As discussed in 
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Sec. | (IE], lattice widths were generally chosen scaled with the square-root of the length. 
Wide systems of W = 16L 1 / 2 were used for most measurements so that W L a for 
any reasonable value of a. This works well for measurements of the current and basin 
fraction which become independent of system width for large widths, but is not suitable for 
measurements of II which approaches unity for any value of F as the system is made very 
wide. 

All of the simulations use toroidal periodic boundary conditions. Any particle leaving the 
left edge of the system appears at the right hand side and any particle leaving the bottom 
row is replaced in the corresponding site in the top row. 



In this paper our primary focus is "sudden" forcing. The lattice is taken to begin with 
all sites below capacity and then the force is suddenly increased to the value of interest. 
This increase causes sites to overflow releasing particles at many sites simultaneously. The 
external force is then held constant and the dynamics of particle flow allowed to proceed. 

This sudden forcing is realized by simply using an (F dependent) initial condition with 
particles placed randomly at each site. Only the number of particles at each site relative to its 
capacity (n) can play a role. We choose this independently for each site, from a distribution 
that, in general, will have weight for both positive (excess) and negative (unsaturated) 
values. Our standard choice is for the initial number of particles at each site [n(x, y,0)] to 
be either +1 or —1 with independent probability p for each site to be above threshold, i.e., 



The parameter p for this particular initial condition is a linear function of the force F used 
in more general discussions. 

If we are studying the steady-state behavior above threshold, then it is sometimes conve- 
nient to only allow sites to start at capacity or above. This removes the site-filling transients 



B. Forcing and initial conditions 




— 1 with probability 1 — p 



1 



with probability p 



(22) 
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and simplifies the behavior; it will be used to study the histograms that show convergence 
to the channel network. 

Although sudden forcing is used for most of the simulations, we will also contrast the 
behavior with that of systems with adiabatic forcing. Adiabatic increase of the force means 
that the force is increased very slowly up to a final value F. Ideally the increase is so 
slow that only one site overflows at a time and any released particle reaches an unsaturated 
site and is trapped before another particle is released. Obviously this is problematic above 
threshold when particles continue to flow indefinitely and interactions between particles must 
be included in some manner. Our simulations with adiabatic forcing will not extend above 
threshold. Below threshold, the absence of collisions between particles and the consequent 
flow of all of them thought primary outlets, enables an efficient algorithm to be used that 



is described in Ref. 20 



III. DYNAMIC CRITICAL BEHAVIOR 

In this section, we present and analyze numerical simulations on the one-deep model 
near to the sudden threshold, focusing on the critical behavior. We make extensive use of 
finite size scaling analysis. 



A. Finite size scaling 

In conventional isotropic systems it is often useful to use an appropriate dimensionless 
quantity that is expected to approach a non-trivial constant value at criticality. For percola- 
tion, the crossing probability of finite size blocks is a convenient choice. A roughly analogous 
quantity in our case is the fraction of systems with a steady-state current I1(L, W, p). In a 
large anisotropic system of length L and width W anisotropic finite size scalingil suggests 
that IT is a function only of the ratios of the correlation lengths to the system dimensions 

n{L,W,p)Mll(L/t f ,W/tfJ, (23) 
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near the critical point. 

If the correlation length in the direction perpendicular to the flow is given by £f± oc 
then we can use an equivalent form 

U(L/^ f ,W/L a ). (24) 

For anisotropic systems it is thus important to know the correct value of a when choosing 
the sizes of the systems if we are to use how II changes with L to obtain useful information. 
Ideally one could perform simulations at a fixed value of 

W 

» - (25) 



Matching finite size scaling expectations to the critical forms of J and B described in Sec. [TC 
for an infinite system, we expect the scaling forms for a system of length L and fixed w to 
be 

Il{L,W,p)nn w (L/£ f ), (26) 
J(L,W,p)*L-V v J w (L/Z f ), (27) 
B(L,W,p)aL- T l" B^L/Zf). (28) 

The w — > oo limits of the current and basin fraction scaling functions should be w- 
independent functions that provide a useful way of making a-independent measurements 
with wide systems. However Yl w (L/^f) approaches 1 as w — > oo for any value of which 
is not useful. 

With the expectation that 

£/ « 7 ~ — \~~ ( 29 ) 

f {p-PcY 



we can write these as 



n(L, W, p) w Ii w ({p - p c )L l / v ) (30) 

J(L, W, p) w L-H v J w ((p- Pc )L^) (31) 

B{L, W, p) w L- T '» P w ((p- Pc )L 1 ^) (32) 
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where 11^, J w and P w are the scaling functions for the particular value of w. As will be 
discussed later, the consideration of the long distance behavior above threshold suggests 
that a should be equal to 1/2 — essentially from random walk scaling of the primary outlet 
paths. We thus first, in the next three subsections, use the value a = 1/2 and the above 
scaling forms to infer the other exponents. Although the fits are good they may be biased 
by the choice of a. 

However, an alternative set of single variable scaling forms can be used without knowing 
a, by working at p = p c . Putting = oo in the form of Eq. (|24]) yields the scaling forms 

U(W,L,p c )^Yl c (W/L a ) (33) 
J(W, L, p c ) « L-^ v J c (W/L a ) (34) 
B(W,L, Pe ) « Lr v l v P c (W/L a ). (35) 

where Il c , J c and P c are single variable scaling functions. These forms will be used, together 
with a measurement of p c , to determine (3/u and Y jv and to try and measure a in Sec. IT I FE. 



B. Scaling of II 

The finite size scaling form of II [Eq. ( [30])] allows us to measure p c and the exponent v 
(at least if we know a) We first identify p c by noting that Il(L,p = p c ) = 11^(0) which should 
be a constant independent of L. This means that curves of YL(p) for different systems sizes 
L should all intersect at p = p c . Fig. f| shows measurements of YL(p, L) for both adiabatic 
and sudden dynamics. Each set of curves intersects quite accurately at a single point. 

The data for II (p) was taken using a series of systems of width W = 4L 1 / 2 . The prefactor 
of 4 was chosen to place the crossing point near II = 1/2 where repeated simulations converge 
most quickly to an estimate of II and where the crossing is most easily seen. By looking 
closer at the crossing of curves for different size systems (as in Fig. [5]) we determine p c = 
0.3133 ± 0.0003 for sudden forcing and p CA = 0.2990 ± 0.0005 for adiabatic forcing. 

With this value of p c we can now examine how the form of II (p) changes with the system 
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size. The finite size scaling hypothesis suggests that a plot of II against L x l v {p — p c ) for 
different lengths should collapse onto a single curve. Such a plot is shown in Fig. [| By 
varying v and comparing how well the data collapses we estimate v = 1.62 ± 0.04. 

Another estimate of v comes from the width of the region over which II (p) is changing.il 
If we differentiate II (p) with respect to p we get a curve which peaks near p = p c and has 
width proportional to L~ x l v . The gradient of a log- log plot of the width of the peak (Ap) 
against L should be —\jv. The width can be obtained from simple integrals of the IT(p) 
data, 



is the center of the peak. No differentiation is necessary if integration by parts is used in 
Eqs. ( |36"D and ([37]), which reduces the noise in the calculations. Fig. ^ shows Ap together 
with a line of the expected slope using the value of v = 1.62 from above. Although the 
data points have the expected linear form they do not give an accurate value for the slope. 
Linear regression gives a rather imprecise estimate, v = 1.69 ± 0.08. 

Our best estimate of v is thus v = 1.62 ± 0.04. This will turn out to be the least 
accurately determined of all our exponents. This is primarily because measuring II is much 
less efficient than measuring B or J. A single run for a large system gives a good estimate of 
B and J, but only provides a zero or one value for IT. Many systems must be averaged over 
to get a good estimate of the fraction that contain moving particles. A further advantage 
of measuring B and J is that wide systems can be used which effects more averaging from 
the larger system and gives results that are less sensitive to the value of a. We thus turn to 
these quantities. 




(36) 



where 




(37) 
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C. Current density 



In an infinite system we expect J ~ {p—PcY- in a finite size system this should only hold 
when £f(p) -C L. Data from simulations of different lengths (and large widths W = 16L 1 / 2 ) 
are shown in Fig. [8] which exhibits the form expected — the data from each system size follows 
a straight line until is of order L. Estimating the slope for an infinite system from this 
figure gives (3 = 1.53 ± 0.03. 

We can also test the finite size scaling form for the current, Eq. (|3"TD, as shown in Fig. |^. 
Varying the value of (3 suggests that (3/v — 0.90 ± 0.03. 

An alternative way to determine (3/v is to use the width scaling form of J [Eq. (|3~4])1- 
This gives a value for (3/v without requiring v to be known (p c is needed but is known more 
accurately). Because the current becomes approximately independent of the width once 
W ^ 2L 1//2 (see Fig. [10]) this gives a rather sensitive test. Working at p c = 0.3133 this gives 



P/v = 0.96 ± 0.01. With v = 1.62 ± 0.04 this translates to (3 = 1.55 ± 0.04. 

Combining this (3/u = 0.96 ± 0.01 value from width scaling with our best value (3 = 
1.53 ± 0.03 gives v = 1.60 ± 0.04 which agrees well with the value of 1.62 ± 0.04 determined 
independently in Sec. pi B . 



D. Drainage basin 



A direct plot of the basin fraction against p — p c (Fig. |TT]) is not as well-behaved as 



the equivalent plot for the current. The range of ordinate values is much smaller and the 
deviation away from the power law is of opposite sign at p — p c small and p — p c large. 
These factors make determining the asymptotic slope rather difficult. The largest sizes are 
approximately linear over one and a half decades of p — p c (compared with almost two and 
a half decades for J).H 

Finite size scaling can also be used for the drainage basin and gives F/u — 0.29 ± 0.02 



(Fig. pi). Varying the width at p c gives Fjv = 0.302 ± 0.005. With v = 1.62 ± 0.04 this 
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translates to T = 0.49 ± 0.02 (see Fig. 



0). 



E. Anisotropic scaling 

As explained above in Sec. |111 A|, it is important to know the value of a in order to choose 



consistent system sizes for simulation. An analogous situation occurs in quantum Monte 
Carlo situations if the dynamic scaling exponent is not known.il In both cases determining 
the anisotropy scaling exponent accurately is quite difficult. 

So far in this section we have assumed that the width of the structures in the system 
scales as the square-root of their length, i.e., that W oc L a with a = 1/2 — although this 
value is not as crucial for the wide, W = 16L 1//2 , systems used for measurements of J and 
B. In this section we provide numerical support for a value close to this. 

Obtaining an accurate estimate of a is rather difficult, a more modest goal will be try 
to rule out a value as different from 1/2 as the value for directed percolation (a ~ 0.63) 

In Sec. [Ill A| we discussed scaling forms for simulations at p = p c which give information 



on the value of a; we should be able to choose the value of a which provides the best data 
collapse with Eqs. fl3"3"|)-(P5|). The problem with this method is that in order to locate p c we 
have to perform a series of simulations using sizes scaling with some value of a. Choosing 
the wrong a will give an apparent value of p c that differs from the value obtained using the 
correct a. This problem is made more acute by the fact that p c is best determined from 
measurements of II on narrow systems which are most sensitive to the value of a. 

To get around this problem we performed additional sets of simulations with different 
system sizes chosen with the values a = 0.4 and 0.6 (the pref actors A in W ~ AL a were 
chosen to give crossings near II = 0.5, but the value of p c is insensitive to the prefactor 



used). As discussed in Sec. |l 1B| we used the crossing of the U w function to give an estimate 
for the apparent critical point for each value of a, p c (ot), obtaining 

p c (a = 0.40) = 0.3150 ± 0.0005, (38) 
p c (a = 0.60) = 0.3110 ± 0.0005. (39) 
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In Sec. |g we found p c (a = 0.5) = 0.3133 ± 0.0003. 

We then performed a series of simulations for different widths and lengths with p chosen 
to equal to the apparent critical parameter p c (a). Plotting II against W/L a should yield 
good data collapse only for the correct value of a [Eq. fl53"D1- These data are shown in Fig. [TJ]. 
The data collapse is most effective for a = 0.5, where it appears to be limited only by the 
statistical error in each point (which is approximately equal to the height of the plotting 
symbols). 

It is also possible to make similar scaling plots for the current and basin fraction. But 
the scaling functions for these quantities rapidly approach a steady-state value as the width 



is increased and so are not very sensitive to the value of a as was seen in Figs. |10| and [13 . 
However, another way to try and fix a is to try to use the other values of p c (a) to fit the data 
for J and B from simulations on wide (W = 16L 1//2 ) samples (which should be insensitive 
to the value of a in a substantial range of a around 1/2) as 16L 1 / 2 > 8L 6 even for L as 
large as 1024). Both the direct logarithmic plots for J and B and the finite size scaling 
plots are significantly worsened if a value for p c as different as 0.311 or 0.315 is used. For 
example, Fig. shows how the logarithmic graph of current is changed by using the a = 0.4 
apparent p c value of 0.315. 

We have provided numerical evidence to support taking a = 1/2. Using similar criteria 
to our other apparent errors, one would guess 

a = 0.5 ±0.05. (40) 

In any case, a directed percolation- like value of a = 0.63 seems unlikely, Again, note that 
most of our measurements were done on very wide systems so the choice of a = 1/2 should 
not be a significant factor if the true a is slightly different. 



F. Adiabatic forcing 



With adiabatic forcing, as mentioned in the Introduction, performing good simulations 
above threshold is too time consuming as the system should be allowed to reach a steady- 
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state before each particle is added. Thus we are limited, essentially, to studies of II, the 
probability that some particles keep moving forever, and the drainage basin fraction at 
criticality. Finite size scaling can be done for Tl(p) for adiabatic forcing (Fig. [TJJ). The 
corresponding estimate for the exponent is v a = 1.60 ± 0.05. If the scaling functions for 
the sudden and adiabatic cases are superimposed from Figs. |B| and [IB] they cannot be dis- 
tinguished. Using the value of the critical point, p c>a , determined from II, Y a jv a can be 
obtained as in Sec. HID. These yield the exponent estimates of Eqs. (|17]) and (plf ) which 



are the same, even within the apparent l-a errors, as those with sudden forcing. 

Note that p c and p CA are rather close; but given the relative smallness of them and the 
observation that with sudden forcing the number of doubly occupied sites after one time 
step is ~ p 3 , it is not surprising that the difference is of this order. Also note that even if 
the fits are done over the range 

\P~Pc\ ~Pc-Pc,a (41) 

to try to separate out possible bias from the critical values being similar, the apparent 
exponents do not change. 

G. Equivalence and scaling 

As mentioned in the Introduction, both our numerical results and analytic arguments 
for the adiabatic forcing suggest a scaling relation between T a and j3 a . Here we explain 



the adiabatic result (following Ref. ^0| for the case of a continuous fluid) and give some 
plausibility arguments for the sudden case which also suggest that both cases are in the 
same universality class. 

If we increase the force adiabatically by an infinitesimal amount 5F then a fraction 5F of 
all the sites will overflow and the current will be increased by the overflow from those sites 
which are on the drainage basin. The current thus increases by 5 J = B 5F, which, when 
integrated, gives the relation 
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Pa = 1 + r a . (42) 

A possible weakness in this argument is that it neglects the anti-correlation between the 
extra sites that overflow and the previous sites that have already overflowed. In the simplest 
version of the model no site will overflow for a second time until all sites have overflowed 
once. This effect could reduce the amount of current induced, but should not be singular 
near threshold and thus should not change the scaling law Eq. fl4"2|). 

If we instead consider sudden increases of the force from an initial value well below 
threshold then the argument is more complicated. If the force is suddenly increased to F, 
there will be a transient current J(F,t) which decays to the steady-state current J(F, oo). 
Now imagine instead increasing the force to a final value of F + SF. The initial transient 
current will now be 

J(F + SF,t = 0) = J(F, t = 0) + 5F. (43) 

The extra particle density will be randomly distributed over the lattice (except for the 
anti-correlation with previously overflowed sites mentioned for the adiabatic case) and some 
fraction of them will survive to join the steady-state current. An extra particle added to any 
site that is on the steady-state drainage basin will increase the final steady-state current by 
one particle, as long as the extra particle passes only through primary outlets (i.e., if it does 
not collide with any other particles), while an extra particle initially off the drainage basin 
will fall into a trap and not contribute to the steady-state current unless it has a collision 
with another particles that forces one of them onto the drainage basin. If we neglect these 
collisions then we would have the same result (/? = 1 + T) as for the adiabatic case with 
the minor difference that the extra steady-state current is not made up entirely of the new 
initial particles: some new particles will fill in traps and allow particles that would otherwise 
have been trapped to survive. 

The complicating factor is that with sudden forcing we cannot neglect collisions caused 
by the new particles. In the adiabatic case the only particles were the new ones with density 
SF so the collision rate was only of order {SF) 2 . In the sudden case collisions between the 
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new particles are similarly negligible but a collision rate of ~ J(F,t)SF between the new 
particles and the existing particles is expected. This is much larger near F c as J(F, t) does 
not vanish for small times. The extra collisions can reduce the steady-state current (which 
will change F c ) but will not induce violation of the scaling relation unless the factor by 
which the steady-state current is reduced is critical. Such a singular contribution is only 
likely to be produced from the long time limit of the decaying current, but in this regime 
the current is small near F c so there will be few collisions and a critical reduction seems 
unlikely. This suggests that (3 = 1 + V holds and also that the presence of splits does not 



change the universality class near threshold. In Sec. [IV D| we will be more quantitative and 
explore the condition on the exponents for this argument to be valid. We will see that it 
also leads to a = 1/2. 

It is possible to relate the exponents we have computed to the fractal dimension of 
the drainage basin on intermediate scales or of large drainage trees below threshold. For 
comparing with other systems we work in general dimensions: 1 downhill and d—1 horizontal. 
Drainage trees of length L smaller than £y contain a typical number of sites ~ L d f . In a 
system of length L ^> £j (with the perpendicular dimensions ~ L a ) we expect the largest 
drainage cluster to contain a number of sites ~ BL 1+a ( d ~ l \ If L <C £/ then the largest 
cluster will contain ~ L d f sites. For these forms to match at L ~ the hyper-scaling 
relation 

d f = l + a{d-i)-^ (44) 



has to hold. In mean field theory, expected^ to be valid for d > 3, it was found that 
df = 4/3, v = 3/2 and T = 1, although the hyperscaling relation [Eq. (f44|)1 will not hold 
except in the critical dimension d = 3 for which logarithmic corrections to mean field theory 
are likely. For d = 2 and taking a = 0.50 ± 0.05, Yjv = 0.28 ± 0.02 Eq. <M) gives 



d f = 1.22 ±0.05. (45) 

There are several useful bounds on the exponents. The generalized Harris criterion for the 
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finite size scaling correlation length for a random system with (d — 1) transverse dimensions 
scaling as L a is, 

v > 2 - (46) 

- 1 + a(d- 1) V ; 

which, with d = 2 and a = 1/2, means v > 4/3. Our measured value of v = 1.62 ± 0.04 
satisfies this bound and is not too far away from saturating it. There is also a simple bound 
on the fractal dimension because the drainage trees must be at least linear, i.e., df > 1, 
implying 

-<a{d-l) (47) 

which is easily satisfied. 

IV. DEVELOPMENT OF THE CHANNEL NETWORK 

The channel network above threshold consists of those sites that carry current in the 
steady-state (see Fig. |I7]). An important feature of our model is that in an infinite system 



the flow converges (albeit slowly) to a fixed network as shown in Ref. [H]. The network is a 
property of the steady-state and is asymptotically time-independent. Which sites are on the 
network depends only on the amount of steady-state current and on the choice of primary 
outlets between sites. The network does not depend on the initial placement of particles, 
how they enter the system or the history of the applied force. 

As discussed in Ref. ^l] these features allow the convergence to the network to be seen 
directly by looking at the number of sites differing between two copies of the same lattice. 
In particular, if particles enter each copy of the system through a different set of sites in 
the top row then we predicted that the number of differing sites should decay as y _1//4 (for 
large y) as the particles move down the system. We have performed simple simulations with 
two such copies and measured a value of —0.262 ± 0.005 for the exponent which agrees well, 
especially considering the difficulties in simulating such a slowly converging quantity. 
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In this section we use a different techniques to follow the development of the steady- 
state channel network via finite size simulations. At the end we will discuss the subtleties 
associated with the slow convergence and the appearance of a second length scale. 



A. Histograms 



In the predicted channel network picture the flow is strictly confined to certain favorable 
channels with other sites containing trapped particles. 

For large systems, the convergence to a fixed channel network means that lattice outlets 
can be divided into two types depending on whether or not they are on the channel network. 
Once a steady-state current pattern has been setup, network outlets pass one particle at 
every time step and off-network outlets never pass particles. In finite systems the network 
is less well defined and the division is not perfect, but we can see the effect developing by 
recording for what fraction of time steps (s) each outlet passes a particle and then plotting 
a histogram for the fraction of outlets with each fractional occupation, h(s). 

As the system size gets large h(s) will approach two delta functions at s = and s = 1 
with the relative weights of each peak determined by the current flowing, 



These peaks are inconvenient to plot so we instead consider the integral of h(s). For an 
infinite system the integral approaches the constant 1 — J/2 for < s < 1. We remove this 
constant by scaling the difference between the integral and 1 by 2/ J and define H(s) by 




(48) 




(49) 



so that H(s = !) = !, 




(50) 



and H should be equal to zero for < s < 1 in the infinite size limit. 
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B. Network scaling regime 



The histogram function H (s, J, L) for large L and small J has a scaling form which is very 
different from that of J, B and II. The correlation length £ n , which is the the characteristic 
vertical scale of the network — roughly the distance between nodes, plays an important role. 
To emphasize this length (and avoid complicating the situation by introducing £f which 
controls the equilibration of the current density) we initially restrict our attention to systems 
that have no traps so that the current is fixed by the initial conditions. The fact that J is 
now an independent variable which we can control is also very helpful. 

A series of simulations were carried out with no traps and initially no more than one 
particle per site. Wide systems of width W = 16L 1 / 2 running for 3L time steps were used 
with histogram data collected for the last L time steps. Each simulation was repeated N s 
times with N s chosen so that WN S = 5120 in order to collect an equal amount of data for 
each size. Basic data for selected J and L values is shown in Fig. [18]. 



With this scaling of the width, the analysis of Ref. [21] suggests the finite size scaling 
hypothesis in the W L 1 / 2 limit, 

Jf(.,J,L) = H.(.,^). (51) 

We expect £ n ~ 1/J 2 for small J so if we plot H(s) against JL 1 / 2 for different values of L 
and J at a single value of s then all the points should fall on a single curve. Plots for three 



values of s are shown in Fig. |19[ Data collapse is quite good with deviations appearing only 
for the largest J values. This shows that the relevant correlation length here appears indeed 
to be £ ra rather than the correlation length seen in finite size scaling in Sec. |iTT| , £/ ~ J~ u ^ . 
We postpone discussion of the role of £/ until later. 

In addition to data collapse for individual values of s we can also consider if as a function 
of s and, with appropriate scaling, collapse all the data onto a single curve (for a limited 
range of J and L). For large systems, L ^> £ n (J), we can expand H n in £ n /L. The fraction 
of time each site spends with the "wrong" occupation should decay as L -1 / 4 (from Ref. PIT) 
so we expect 
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H(s,L^UJ))*Hl(s) (|) ' . (52) 

since in the limit L — > oo we expect = 0. This means that LV*H(s, J, L) should 

collapse onto a set of curves for different J as L — > oo. This works rather well as seen in 
Fig. ^0] which shows scaled results for 16 different pairs of J and L values which overlap 
onto five different curves for the five different values of J. Overlap is least effective for 
the smallest J value which corresponds to the largest £ n (J) where the limit L ^> £ n is not 
reached. 

Finally we also know that £ n (J) ~ J~ 2 for small J, so H(s,J,L) should be the same 
for the same values of J 2 L for J<C 1. This works quite well, and when combined with the 
L-scaling we expect it to produce collapse onto a single curve in a limited regime: 

( J 2 L) 1/4 #(s, J, L) « H n (s) if J < 1 and J 2 L > 1 (53) 

The collapse is quite good, a plot is shown in Fig. ^lj, although the regime of applicability 
is rather small as we also need the number of particles in the system ~ JL 3//2 to be not too 
small in order to have sufficient data. 

C. Critical scaling regime 

So far, we have only investigated the local current distribution in the regime where 
L ~ £ n ~ 1/J 2 or larger, i.e., the regime important for formation of the steady-state 
network. But from the data of Fig. ^ we can see that the mean current density has converged 
to its large system — equivalent to long time — limit when L w £j £ n . What then is the 
distribution of local current densities in this critical scaling regime? 

A simple scaling argument suggests the answer: In a wide system of length L ~ £j (with 
periodic boundary conditions), each section of width £f± ~ will have of order one (but 
sometimes zero as can be seen from the behavior of IT) current paths in the steady-state. 
Thus a fraction of order of the sites will carry a steady-state current. Since the total 
current density is J ~ £/ j the typical (time averaged) current through theses sites will be 
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8 ~ ji-W/3. (54) 
We can thus guess a scaling form for the scaled cumulative distribution H(s) in this regime: 

H & 1 L ) - JI^Hf (71^ LJV ' P ) ( 55 ) 

where for simplicity we have again simplified to the large W (3> L a ) limit. For convenience 
we have used J instead of p — p c as the scaling variable, although, of course now J is 
determined by the dynamics — including collisions and trap-filling — on length scales ^ £/. 



D. Intermediate regime and scaling form 

For consistency we expect that for large values of the scaling arguments Hf should match 
with the network scaling function H n [Eq. (|5lD l in the limit of small values of the latter's 
scaling arguments. Thus for 

^ < L < fn (56) 

we expect a single argument scaling function whose argument must be a product of the 
scaling arguments in both regimes. This yields a unique choice of the argument: 

H(.,J,L)«\H.yL ? ) (57) 

with 

In the intermediate regime, neither collisions nor critical effects should play much role. 
Thus we should be able to understand the scaling of Eq. ( [5T|) from simple considerations. 
Regions of length i ^> £y will contribute current to the drainage basin and hence to the 
steady-state current unless they happen to be the exponentially rare regions which sit in 
anomalously large holes in the drainage basin. Furthermore, regions much further apart 
than ^ (or will contribute roughly independent currents. Thus on scales 3> £/, the 
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current will collect on drainage trees with contributions from each region of length i 3> £/ 
being simply J times its area (with small variations around this). As soon as one particle 
can pass down a drainage tree, then this branch must have no traps so that traps cannot 
play a role in the behavior of large drainage trees. Also, by assumption — to be verified 
later — collisions are not important in this regime. The statistics of the drainage trees, the 
current network, and the local current density on them in this intermediate regime can thus 
only depend on J, L, and properties of the primary outlet trees; they are therefore entirely 
determined by random walk statistics. The expected scaling for local current densities is then 
very simple: the periodic boundary conditions from bottom to top imply that the current 
network consists of the primary outlet paths from all points in the top row which emerge 
at the same point in the bottom row. These will be separated by horizontal distances of 
order L 1//2 and thus have basins of width L 1//2 implying time averaged local currents on the 
network of typically 



s ~ 



JL 1 ' 2 (59) 



with a distribution determined by the random walk properties. This is consistent with 
Eq. ( ]57| ) only if e = 1. This implies that a = 1/2 exactly. Indeed from the above discussion, 
we could have guessed this from the observation that the system in this regime does not 
"care" about the physicslil'il that determines /3/u. 

To check that the above argument is consistent, we need only to check that collisions are 
rare in this intermediate regime: since 

s ~ JL 1 ' 2 < 1, (60) 

this immediately follows. 

Results from simulations that try to reach this intermediate regime are shown in 
Figs. ^ and|23|. Initial conditions that include traps (fixed p) and exclude traps (fixed J) are 
compared for similar mean J. Fig. ^2] shows H(s) for the two cases averaged over 28 systems 
of a single size. The parameters chosen correspond to J 2 L « 0.1 and ~ 2.5 [estimating 
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£/ ~ (p — p c ) v with J m 2(p — p c Y from Fig. §]. To be in the intermediate regime we need 
J 2 L -C 1 and L £f, if we neglect the difference between v and (3 then this means we need 

J 2 L < 1 « JL. (61) 

This is difficult to achieve because the minimum accessible current density is limited by 
the size of the system. Despite this difficulty, the data collapses to the scaling form of 
Eq. (f>7\ ) quite well (Fig. ^3[) . However, the two different initial conditions do not appear 
to have converged to the same scaling function in the accessible range of the parameters. 
Nevertheless, our conjectures on the intermediate regime would seem to be consistent with 
the numerics. 

The existence of an intermediate regime with L ^> £f but collisions still unimportant, is 
closely linked to the arguments of Sec. |III G| that collisions are unimportant at long times 



near threshold. To see this, consider a system just below threshold with 5 = F c — F and 
correlation length £/ ~ 5~ u , which is the length of the longest typical drainage tree. If the 
force is suddenly increased by, say, 8/2, a number of particles of order 8£/ will overflow on 
a tree of length £ j resulting in an average number of particles appearing at the same time 
at the bottom of the tree of order 

n s ~ 8£ d /~ X . (62) 

(They will arrive over a time interval £/.) The natural guess is that the condition for 
collisions to be irrelevant at the critical point is that for £j large, we must have n s <l, i.e., 
that 

1 - v{d f - 1) > 0. (63) 



Using the equalities of Sec. Ill (j[ df = l + a — T/v, a = 1/2 and /3 = 1 + T, this is equivalent 
to 
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which is identical to the condition that £ n ^> £/ just above threshold. We thus conclude 
that because of the inequality [Eq (0)] collisions are dangerously irrelevant at the critical 
point — but crucial for the flowing phase at long times; that adiabatic and sudden forcing are 
in the same universality class; and that all the scaling laws should be correct leaving just 
two independent exponents {(3 and v). In three dimensions, where mean field theory results 
should hold up to logarithms,^ one similarly finds that collisions should be dangerously 
irrelevant near but above threshold. 



V. APPLICATIONS AND CONCLUSIONS 

This paper has investigated a new universality class of dynamic critical phenomena, 
which appears to be rather robust. Perhaps the most interesting aspect of the scaling is 
the presence of two lengths above threshold, or, equivalently, two different time scales for 
equilibration of the current and of the flow pattern (£ n ). We showed how the histograms 
of current distributions indicated the presence of the longer length (£ n ) and how to reconcile 
this scaling (with £„,) with that of £/. The other surprising result is that the critical behavior 
seems to be in the same universality class regardless of whether the force increase is sudden 
or adiabatic in spite of the dependence of the critical force on the history. This observation 
and analytic arguments suggest that dynamic particle collisions are dangerously irrelevant: 
they do not effect the critical exponents unless they are completely absent in which case the 
behavior is very different with all of the particles becoming concentrated onto a single large 
tree at very long times. 

The model we have studied updates all of the sites synchronously so that particles on 
different rows can never interact or intermingle. A more realistic model would allow some 
inter-row diffusion. The fact that sudden and adiabatic forcing seem to give the same 
critical behavior suggests that the exponents are also likely to be independent of the order 
in which the particles are moved and other time-independent changes in the local rules such 
as allowing for position or occupation dependent local particle velocities. However allowing 
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inter-row diffusion will change the long time approach to the network (on times greather 
than £ n ) from t -1 / 4 to t -1 / 2 as explained in Ref. |21]. 

Previously Narayan and FisherS studied a related model with a continuous fluid instead 
of discrete particles and adiabatic force increase, obtaining numerical results below threshold 
and using scaling arguments to infer the properties of the current network above threshold. 
They measured 

d f = 1.21 ±0.02 (65) 

using scaling and assuming that a = 1/2. This is very close to our value for discrete particles 
(df = 1.22 ± 0.05, the larger error arising from the uncertainty in a, if a = 1/2 is assumed 
exact then our error is also ±0.02). They also obtained 

v= 1.76 ±0.02 (66) 

which can be compared to our value of v — 1.62 ± 0.04. The difference is twice the apparent 
errors but, as the measurements were made by different methods, agreement is certainly 
reasonable. We conclude that the below threshold behavior (and thus £/) of both continuous 
and discrete particles are in the same universality class, 

But above threshold the models differ because the continuous model has a different 
behavior of £ n , since the average depth of the rivers approaches zero at threshold. In the 
continuous fluid model, £ n depends on the behavior of the probability density, p(b), that 
the secondary outlet barrier is a small amount, b, higher than the primary barrier. For 
p(b) ~ const, as b — > 0, £ n <C near threshold. This implies that, in contrast to our 
discrete model, the river splits must play a role near threshold if the force is increased 
rapidly, a situation not considered in Ref. |20|. 

We conclude with a few comments on possible connections between the critical phenom- 
ena discussed here and vortex motion in dirty superconducting films. 

One possible behavior of vortices that is very different from that found in our model is 
for the time averaged vortex current in most regions of space to be non-zero even near to the 
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critical current threshold. Because of strong tendencies of a vortex lattice to break up near 
threshold, this motion would have to involve sections of vortices moving together but not 
synchronously with neighboring regions. It may be possible that such non-uniform irregular 
motion with different regions moving at different times could persist in steady-state. But 
this could also be a transient phenomenon with the system eventually settling down, near 
the critical current, to a steady-state pattern of channels separated by wide regions with 
vortices at rest. If this is the case, it is quite possible that, given the degree of robustness 
of the critical behavior found here so far, the critical current phenomena would be in the 
same universality class as our model even with the complications of longer range vortex 
interactions, particles stopping and starting, history dependence of the critical current, etc. 

A signature of the channel behavior could be found by constructing a histogram of 
the time averaged local vortex velocities which should show a large peak at zero velocity 
near threshold and some distribution — with small total weight — at non-zero velocities, in 
the simplest scenario centered on a velocity that does not vanish at threshold. Behavior 
qualitatively like this has been seen in simulations of vortices in dirty thin films.0 Other 
possible experimental probes of channel flow will be discussed elsewhere. 

It should be noted that the critical behavior found here is restricted — at least in the 
model — to the first increase of the force. Decreasing the force and subsequent increases can 
result in different behavior. This and other qualitative phenomena that can occur in these 
types of models and vortex systems will also be discussed elsewhere. 
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FIGURES 
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FIG. 2. Fraction of systems with a non-zero steady-steady current [n(p)], steady-state current 
[J(p)] and fraction of sites in the drainage basin [B(p)] for a 1024x512 lattice averaged over 32 
systems. 




FIG. 3. Lattice model below threshold. The black sites are saturated and the white sites are 
unsaturated. The lines show the primary outlets of each saturated site and form drainage trees. 
The force direction is down the page. 
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FIG. 4. Fraction of systems that have a steady-state current for sudden and adiabatic increase 
of the force. Data use the same set of system sizes with widths W = 4L 1 / 2 and using 512 samples 
for each. 
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FIG. 5. Close up of II(p) curves for selected lengths with sudden forcing. Error bars are from 
measured sample variance. Crossing point is estimated as p c = 0.3133 ± 0.0003. 
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FIG. 6. II finite size scaling with p c = 0.3133 and v = 1.62 for systems of width W = 4L 1//2 
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FIG. 7. Width (Ap) of peak in dU/dp for systems of different lengths, and widths W = 4L 1 / 2 . 
Line has slope —\jv with v = 1.62. 
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FIG. 8. Log-log plot of J versus (p — p c ) for p > p c with p c = 0.3133. Line shown has 
slope P = 1.53 
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FIG. 9. Finite size scaling plot for J using p c = 0.3133, v = 1.62 and (5 = 1.53. 
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FIG. 10. Scaling of current at p = 0.3133 with sample size, using (3/v = 0.96 




FIG. 11. Log-log plot for basin density. 
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FIG. 12. Finite size scaling plot using p c = 0.3133, V jv = 0.29 for the basin fraction B. 
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FIG. 13. Scaling of drainage basin fraction at p = 0.3133 with sample size using T/u = 0.302 
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FIG. 14. Tests for the width scaling for different values of a. Simulations run at p c (ct). Data 
for a = 0.6 is offset vertically by +0.2, data for a = 0.5 is offset by +0.1. 
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FIG. 15. Log- log plot of current using the value p = p c (ce) determined from using widths scaled 
with a = 0.4. Compare with Fig. S 
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FIG. 16. II finite size scaling for adiabatic increase with p Cja = 0.299 and v a = 1.60 and using 
narrow (W = 4L 1 / 2 ) lattices 
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FIG. 17. Channel network showing which outlets in the lattice carry current in the steady-state 
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FIG. 18. Basic cumulative histogram data from one-deep model using J = 1/16, 1/4 and 1 
with system sizes L = 128, 512, 2048. For a given L, the curves for larger J are flatter. 
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FIG. 19. Scaling of histograms for fixed occupancy fraction (s) with system size (L) and current 
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FIG. 20. Length scaled histogram data showing overlap for different values of J for the large L 
limit. Currents used are J = 1, J = 1/2, J = 1/4, J = 1/8 (only for L > 256) and J = 1/16 (only 
for L > 512) 
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FIG. 21. Data collapse for selected parameters: J = 1/8 for L = 256,512,1024,2048 and 
J = 1/16 for L = 1024,2048. 
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FIG. 22. Histograms for current distribution from initial conditions with and without traps. 
Using a system of size 512x360. System with traps has p=0.35 and an average current of J = 0.0134. 
Current without traps is chosen to be J = 0.0134. 
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FIG. 23. Test of the scaling form for H(s) in the intermediate regime [Eq. (| 
of system sizes from L = 128 to L = 1024 with p and J chosen so that L 2 
0.02 < s < 0.10. 
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